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ABSTRACT 

Context. In the framework of the coagulation scenario, kilometre-sized planetesimals form by subsequent collisions of pre- 
planetesimals of sizes from centimetre to hundreds of metres. Pre-planetesimals are fluffy, porous dust aggregates, which are in- 
homogeneous owing to their collisional history. Planetesimal growth can be prevented by catastrophic disruption in pre-planetesimal 
collisions above the destruction velocity threshold. 

Aims. We assess whether the inhomogeneity created by subsequent collisions has a significant influence on the stability of pre- 
planetesimal material to withstand catastrophic disruption. We wish develop a model to explicitly resolve inhomogeneity. The input 
parameters of this model must be easily accessible from laboratory measurements. 

Methods. We develop an inhomogeneity model based on the density distribution of dust aggregates, which is assumed to be a Gaussian 
distribution with a well-defined standard deviation. As a second input parameter, we consider the typical size of an inhomogeneous 
clump. For the simulation of the dust aggregates, we utilise a smoothed particle hydrodynamics (SPH) code with extensions for mod- 
elling porous solid bodies. The porosity model was previously calibrated for the simulation of Si02 dust, which commonly serves as 
an analogue for pre-planetesimal material. The inhomogeneity is imposed as an initial condition on the SPH particle distribution. We 
carry out collisions of centimetre-sized dust aggregates of intermediate porosity. We vary the standard deviation of the inhomogeneous 
distribution at fixed typical clump size. The collision outcome is categorised according to the four-population model. 
Results. We show that inhomogeneous pre-planetesimals are more prone to destruction than homogeneous aggregates. Even slight 
inhomogeneities can lower the threshold for catastrophic disruption. For a fixed collision velocity, the sizes of the fragments decrease 
with increasing inhomogeneity. 

Conclusions. Pre-planetesimals with an active collisional history tend to be weaker. This is a possible obstacle to collisional growth 
and needs to be taken into account in future studies of the coagulation scenario. 

Key words, hydrodynamics - methods: numerical - planets and satellites: formation - planets and satellites: dynamical evolution and 
stability - protoplanetary disks 



1. Introduction 

Terrestrial planets are thought t o form by core accretion in pro- 
toplanetary accretion discs (e.g. lLissauer & StewartI [19931) . The 
discs consist of both hydrogen and helium gas and a small frac- 
tion of dust , which in the early stag es exists as micron sized dust 
grains fe.g. lDuUemond et al.ll2007f) . Interaction between the gas 
and dust induces relative motions bet ween the dust grains and 
causes collisions among them (e.g. IWeidenschilling & Cuzzil 
Il993h . Mutual sticking by van der Waals forces occurs and 
by this process millimetre- t o centimetre-sized, h ighly porous 
pre-planetesimals form (e.g. iBlum & WurmI l2008h . Additional 
growth by means of collision is endangered by several obstacles: 
rebound, mate rial loss by accretion or photoevaporation, and 
fragmentation. IZsom et al.l (1201 Ol) found that the growth of pre- 
planetesimals to kilometre-sized planetesimals may be halted 
by mutual compaction and bouncing, which leads finally to a 
populatio n of dust and pebble-siz ed objects in the protoplane- 
tary disc. IWeidenschillind (Il977h found that metre-sized bod- 
ies have the highest drift speed towards the host star, where the 
objects get accreted or photoevaporated. For these objects, the 
drift timescale is much shorter than the time needed for a planet 
to form. One of the most serious obstacles to planetesimal for- 
mation is the fragmentation barrier: with increasing aggregate 
size the mutual collision velocities also increase and catastrophic 



disruption might become more common than mass gain (e.g. 
iBrauer et al. 2008). 

If by some mechanism a sufficient population of kilometre- 
sized objects is formed, their self-gravity acts as an accretion 
mechanism and ensures further growth to planets. In this gravity 
dominated regime, kilometre-sized objects and yet larger ones 
are referred to as planetesimals. 

To investigate p lanetesimal formation, global dust coagu- 
lation model s le.g. iDuIlemond & DominikI I2005L IZsom et aH 
,2010, Birnst iel et al .1120101) numerically compute the evolution 
of dust aggregates in a protoplanetary disc starting from micron- 
sized dust grains. These models require detailed collision statis- 
tics evaluated for the outcome of two-body pre-planetesimal col- 
lisions as functions of important parameters such as mass ra- 
tio, porosity, impact parameter, and relative velocity of the col- 
lision partners. In the past decade, many expe rimental studies 
have bee n carried out to collect these data (e.g. iBlum & WurmI 
2008, Gu ttler et alJl20ld for summaries). However, laboratory 
setups become difficult and infeasible for centimetre-sized pre- 
planetesimals and larger ones since they have to be carried out 
under protoplanetary disc conditions, i.e. in a vacuum and under 
micro-gravity. In this regime, collision data has to be acquired 
by numerical simulations. 

In their pioneering work, iBenz & Asphaugj 11991 [T99I de- 
veloped a solid body smoothed particle hydrodynamics (SPH) 
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implementati on for brittle material based on the dynamical frac- 
ture model bv lGradv&Kippl lllQStf). which evolves the propaga- 
tion of cracks starting from an init ial Weibull ( 1939) distribution 
of flaws. On these grounds, lBenzl(l2000l) investigated collisions 
between rocky non-porous pre-planetesimals in their typical col- 
lision velocity regime from 5 ms ' to 40ms~' and found o nly 
disruption. However, experiments (e.g. lBlum & Wurmll2008h in- 
dicate that pre-planetesimals are fluffy, porous objects rather 
than solid rocks. For this reason. ISironol ( |2004|) developed an 
SPH implementation of a porosity model based on porosity- 
dependent material quantities such as bulk and shear mod- 
uli as well as compressive, shear, and tensile strengths. In his 
simulations of porous ice aggregates, he found that owing to 
their porous structure the kinetic energy can be dissipated ef- 
fectively, permitting collisional growth. Sirono's approach in- 
cluded a dama ge model base d again on the dynamical frac- 



ture model of Gradv & KippI (Il980h . In addition, a damage 
restoration mo del tries to capture d amage healing by compres- 
sion. However. ISchafer et al.l (l2007h found that Sirono's damage 
model is not applicable to their SPH simulations of porous ice 
because it includes damage by compression. In contrast, mate- 
rials such as porous ice and porous dust form new molecular 
bondings, when they are compressed. Thus, their strengths in- 
crease but do not decrease as predicted by t he Sirono damage 
model . On the basis of the implementation of 'Benz & Asphaud 
(11994 [1995) . Jutzi et al. (2008) combined the damage model of 
iGradv & KippI ( Il980l) and the p-a model of E errman ril ([T969h 
to simulate porous brittle media. The model was calibrated wit h 
laboratory high velocity impact experimen ts dJutzi et al.ll2009h . 
However, the model by lJutzi et all (|2008|) is also not applica- 
ble to porous protoplanetary dust because it includes damage by 
compression and has no suitable model of damage healing. 

To simulate SiOa dust, which is used as protoplanetary dust 
analogue in many experimental stud ies, we develo ped a suitable 
porosity model based on the ideas of lSironol (l2004h . that was im- 
plemented into a solid body SPH code to numerically investigate 
pre-planetesimal collisions dGeretshauser et alll2010|). With the 
aid of laboratory benchmark experiments (see also Giittler et al.l 
l2009l) . the porosity model was calibrated and tested for the sim- 
ulation of Si02 dust. In particular, the simulation reproduced the 
compaction, bouncing, and fragmentation behaviour of homo- 
geneous dust aggregates of various porosities in a quantitatively 
correct way. However, this approach does not include an explicit 
damage model. 

In this paper, we present an approach to explicitly resolve in- 
homogeneity of dust aggregates. Since an inhomogeneous den- 
sity distribution is one of the main sources of fracture in porous 
aggregates, this work mimics a damage evolution by the evo- 
lution of the density. The inhomogeneity is described by a 
Gaussian density distribution and the typical size of a clump. The 
aggregate inhomogeneity of SiOa dust can easily be measured 



in the laboratory by X-ra y tomography ( see Guttler et al .1120091 
Fig. 10). In contrast, the loradv & KippI (1 19801) damage model 
relies on the correct parameters for the initial Weibull distribu- 
tion of flaws, which are difficult to determine in laboratory ex- 
periments. Our approach of resolving inhomogeneity is closely 
connected to the porosity model, where the strength quantities 
are porosity dependent. Therefore, any fluctuation in porosity 
induces a fluctuation in strength over the aggregate, which is 
utilised to model the evolution of the material under stress. 

To classify the collisional outcome of our simulations 
involving inhomogeneous aggregates, we utilise the four- 
population model. This cl assification scheme was introduced by 
iGeretshauser et al](l201 ll) to enable a precise mapping of all col- 



lisional outcomes involving any combination of sticking, bounc- 
ing, and fragmentation processes. Four fragment populations are 
distinguished according to their masses: both the largest and the 
second largest fragment as well as a set of fragments whose mass 
distribution can be approximated by a power law. To take into ac- 
count the resolution limit, the sub-resolution population contains 
all fragments that consist of a single SPH particle. 

The outline of this paper is the following. In Sect. |2] we 
describe the SPH method with extensions for the simulation of 
solid bodies and the applied porosity model. Section|3]firstly de- 
scribes a previous study of a damage model, then we present our 
approach based on the inhomogeneity of dust aggregates and its 
implementation. The results of collisions with inhomogeneous 
aggregates are discussed in Sect. |4| We summarise our findings 
in Sect. |5] where we also discuss the results and provide an out- 
look on future work. 



2. Solid body SPH and porosity model 

For our simulations, we apply the smoothed particle hy- 
dro dynamics (S PH) numerical method originally developed 
by iLucyl ( 1 19771) and |Gingold & Monaghan (1977) for astro- 
physical flows. In subsequent studies, SPH was extended to 
model th e elastic and plastic deformations of solid materi- 
als (e.g. iLiberskv & Petschekl 1199 iL iBenz & Asphaugj 11994 
iRandles & Liberskv' '\996). Our parallel SPH code parasph 
tHipp & Rosenstie l 2004, Schiifer 2005, Schafer et al. 2003 
IGeretshauser et all[2010l) is based on these concepts. In the SPH 
scheme, the time evolution of a solid body is simulated by means 
of the Lagrangian equations of continuum dynamics 
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which represent the continuum and momentum equations, re- 
spectively. The quantities p, i',,, and cr„^ have their usual mean- 
ings, i.e. density, velocity, and the stress tensor The Greek in- 
dices denote the spatial components and the Einstein sum con- 
vention applies to them. Since our equation of state (EOS) is 
energy-independent, we do not solve the energy equation, hence 
omit to state it here. Within the framework of SPH, the con- 
tinuous quantities are discretised into mass packages, which we 
refer to as "particles". These represent the sampling points of 
the method and interact with each other within a limited spatial 
range, called the smoothing length h. The influence of a particle 
^7 on a particle a depends on the particle distance |x" - x''| and is 
weighted by the kernel function W(|x" - 'sl'\,h). This is usually 
a spherically symmetric function with compact support, which 
is differe ntiable at least to first order. W e utilise the cubic spline 
kernel of iMonaghan & Lattanzid (Il985l) . 

The particles of the SPH scheme move according to the 
Lagrangian form of the equations of continuum mechanics. 
Thus, they represent a natural frame of reference for any defor- 
mation of the simulated solid body. In the context of fragmen- 
tation, two fragments are considered to be separate once their 
constituent subsets of particles no longer interact, i.e., when the 
fragments are separated by more than 2 h. The fragments may be 
spread out over an unlimited computational domain. This is why 
computations are carried out at the particle positions. 

The stress tensor in the momentum equation can be divided 
into the term accounting for pure hydrostatic pressure p and the 
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term representing pure shear, which is given by the traceless de- 
viatoric stress tensor 5 Hence, 



(3) 



In contrast to viscous fluids, p and Sap do not depend on the 
velocity gradient but on the deformation of a soUd body given 
by the strain tensor 



2 \ dxP dx" ) ' 



(4) 



where x' are the coordinates of the deformed body. In particular, 
the relations that hold are 



5"^ = 2//Je"^- ^(^"V'J- 



rotation terms, 



(5) 
(6) 



where the quantity d denotes the dimension and K and ji are the 
bulk modulus and shear modulus, respectively, for elastic de- 
formation. To obey the principle of frame invariance, rotation 
terms have to be added to the deviatoric stress tensor. For this, 
we apply the Jaumann rate for m, which was previo usly used 
by iBenz & Asph"aug ('1994|) a nd ISchafer et alJ (l2007h . and de- 
scribed bv .Geretshauser et all (l201d) in more detail. 

The above relations (Eqs.|5]and|6ll describe the time evolu- 
tion of a solid body in the elastic regime, i.e., for reversible de- 
formations. With respect to pre-planetesimal collisions, we are 
concerned with highly porous material that undergoes elastic and 
plastic deformatio ns. Both aspects are described by means of the 
porosity m odel of iGeretshauser et all (|2010|) . which is based on 
the ideas o f lSironol (120041) ' and can be categorised as a plasticity- 
based porosity model. 

In this approach, the quantities for the elastic regime (bulk 
and shear moduli), as well as the strength quantities (compres- 
sive, shear, and tensile strength), depend on the filling factor, 
which accounts for the sub-resolution porosity. The filling factor 
<p is related to the density p and the porosity (D by 



Ps 



(7) 



where ps is the density of the matrix material. We use a con- 
stant value of Ps - 20 00 kg m~^ for our SiOi material (e.g. 
iBlum & Schraplej|2004l) . In particular, the elastic quantities are 
related to the filling factor by a power law 



K{(l,) = 2^i((|,) = 



(8) 



where y - A,Kq is the bulk modulus of an uncompressed random 
ballistic deposition (RBD ) dust sample with 0rbd = 0.15 (e.g. 
iBIum & Schrapleill2004t) . such that Kq = /:(<^rbd) = 4.5 kPa. 

For the compressive, tensile, and shear streng t h we 
adopt some empirical rel ations (see 'Giittler et alT l2009l 



aaopt some empirical rel ations (see uuttier et aiJ izuu^ 
IGeretshauser etalllloTol and IGeretshauser et al.i .201 ll for de- 
tails). These relations allow us to reproduce the compaction, 
bouncing, and fragmentation behaviour to high accuracy in a ve- 
locity range of the order of 0.1 to 10 ms ' (Geretshaus er et alj 
I2OIOI) . IGeretshauser et alJ (1201 ll) showed that all collision types 
found in laboratory experiments can be reproduced. In the same 
reference, the laboratory velocity thresholds for the transitions 
between sticking, bouncing, and fragmentation behaviour re- 
sembled the simulation results for collisions of medium porosity 
aggregates with up to 27.5 ms"'. Therefore, we expect these re- 
lations to be valid for collisions below the sound speed of the 



dust material, which is ~ 30 ms 1 jGiittler et al.ll2009h . For the 
compressive strength, 

X A In 10 



Z(0) = p„ 



^max V^min 



0max - <A 



- 1 



(9) 



where 0n,in + s < (p 

^min -0.12 and ^^nax 



< 0max and s = 0.005. The quantities 
- 0.58 are the minimum and maximum 
filling factors, respectively, in the compressive strength relation. 
Both values can be exceeded by the material. The power of the 
compressive strength relation is In(lO) times the parameter A 
with A = 0.58. The quantity p^ = 260 Pa accounts for the 
mean pressure of the compressiv e strength relation. According to 
'Geretsha user et all (1201 Oi 1201 ll) . these parameters are valid for 
im pacts below ~ 30 nis~' . Th e influence of A and p^ was studied 
bv IGeretshauser et al.l (l2010l) . For (p < ^min + s, the compressive 
strength relation is continuously extended by the constant func- 
tion Z(0) = £(0min + e) and for i^niax ^ <P, we set E(0) - 00. The 
tensile strength is given by 

r((^) = -lO^-^+'-'^^'^Pa (10) 

and the shear strength is the geometric mean of the compressive 
and tensile strength 



(11) 



We note that the three strength quantities are filling-factor de- 
pendent. This is because as the filling factor increases the 
monomers of the dust material are packed more closely. Thus, 
each monomer establishes more bonds, which increases the 
strength. The filling-factor dependence of the strength quantities 
is exploited by the approach presented in this paper Together the 
bulk modulus, the compressive strength, and the shear strength 
yield the full equation of state for pure hydrostatic compression 
or tension 



f2(0) 



1) 



for (p^ < 
for <(p <<p^ 
for (p < 0J 



(12) 



where (p^ > (p^, and (p^ and 0J are critical filling factors. The 
value of (p^ marks the transition between elastic and plastic com- 
pression and (p^ defines the transition between elastic and plastic 
tension. The filling factors in-between the critical values repre- 
sent the elastic regime. The quantity 0^ is the reference filling 
factor that represents the filling factor of the porous material at 
vanishing external pressure. Once the material undergoes plastic 
deformation, the reference filling factor is reset using the mid- 
dle relation of Eq.[T2l Details of th i s pres sure reduction process 
can be found in IGeretshauser etal] ( l2010l) . In the elastic regime 
(middle line of Eq. \T2\ . the pressure varies with density accord- 
ing to the elastic path given by Eq. [8] If 0^ is exceeded, the 
pressure follows the much shallower compressive strength path 
2(0) in the positive 0-direction, whereas if (p^ is exceeded, the 
pressure follows the shallower ten sile strength path T((p) in the 
negative ^-direction (cf. Fig. 1 in IGeretshauser et"al]|2010l) . In 
both cases, the absolute value of the pressure, which is com- 
puted according to the elastic equation of state, is reduced to the 
compressive or tensile strength, respectively. This accounts for 
the stress relieved by plast ic deformation. Details of the imple- 
mentation are described bv IGeretshauser et alj ( l2010h . 

For the plastic deformation by pure shear, the deviatoric 
stress tensor S „fi is reduced according to the von Mises plasticity 
given by (e.g. lBenz & AsphaudI 19941) 



(13) 
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The function / is defined by / - min |y^/3/2, l]- A n alternative 
approach (e.g. IWilkinsl[l964 iLiberskv et al.|[T993l) uses the ex- 
pression / = min 1^ -\/Y^/3J2, l] but we expect this to lead to dif- 
ferences from our adopted approach that are marginal. In this re- 
lation, J2 - 0.5 S apS ap and Y - Y{(p) are the second irreducible 
invariant of the deviatoric stress tensor and the shear strength, 
respectively. 

3. Resolving inhomogeneity 

3.1. Previous approach 

In contrast to ductile media, brittle materials, such as basalt, 
granite, or porous pumice, do not rupture by plastic flow. This 
is because the material is not completely homogeneous but con- 
tains little flaws. These are little defects in the medium. With in- 
creasing strain, cracks develop originating from these flaws and 
start to pervade the solid body. In brittle media, stress is relieved 

by developing cracks. 

To model damage, iGradv&KippI (119801) introduce a scalar 
parameter Q < D < 1, where D = represents undamaged and 
D = 1 fully disintegrated material. The stress tensor (Eq. O is 
modified according to 



O-dam = (1 -D)o- . 



Thus, damaged material may feel less stress than undamaged 
material. The local damage D is defined as the ratio of the vol- 
ume defined by the growing crack to the volume in which the 
crack is growing (see lGradv &^iDDllT980h . The local time evo- 
lution of the damage D is given by 



dt 



(15) 



where Wact accounts for crack accumulation and denotes the num- 
ber of activated flaws, and is the radius of a circumscribing 
sphere in which the crack evolves. It accounts for the maximum 
size of a damaged area. The quantity Cg is the crack growth ve- 
locity. The concept of flaw activation originates from the idea 
that cracks do not start to grow for any strain applied but are ac- 
tivated for some strain threshold eact- The number of flaws n per 
unit volume with tact is given mostly by a power law 



"(fact) 



«^wbtact 



(16) 



This dis t ributio n of flaws in a brittle material was proposed by 
IWeibullI (Il939h . It is based on two material parameters fc„b and 
niv/b, where k„b is the number of flaws per unit volume. The flaw 
distribution according to Eq.[T6]is set as an initial condition for 
the material before the simulation. 

A disadvantage is that the material parameters ^^b and 
mwb are rarely available because they are difficult to measure. 
Unfortunately, small variations in A:„b and ;7Zwb lead to large dif- 
ferences injheacth^at^ and simulation results (e.g. 
ISchafeill2005llJutzi et al.ll2009h . 

Because the damage model by iGradv & KippI (Il980l) is 
only suitable for brittle material, it is not applicable to pre- 
planetesimal material represented by Si02 dust, which has prop- 
erties that are between those of a ductile and a brittle material. 
Therefore, it is desirable to develop a simple damage approxima- 
tion consistent with the applied porosity model. To summarise, 
three ingredients are essential for a damage model: (1) a distri- 
bution of defects, (2) the stress reduction due to damage, and (3) 
the typical size of a damaged area. 



3.2. The approach of resolving inhomogeneity 

In contrast to the damage model presented in Sec. 13.11 our ap- 
proach does not explicitly consider the evolution of defects, 
hence is not strictly a damage model. However, it is designed to 
behave as an approximation to a damage model in some aspects, 
which are described in this section. As a significant advantage, 
our approach relies on quantities that that can be more easily 
measured by laboratory experiments than the parameters of the 
model by Grady & Kipp ( 1980). 

The starting point for our approach is the inhomogeneous 
nature of the material. Macroscopic dust aggregat es created by 
the R BD method are not completely homogeneous dGiittler et al.l 
'2009'). The filling factor instead was found to follow a Gaussian 
distribution around a median of 0^ ~ 0.15 with (p^- = 0.013. 

According to the porosity model of Sec. |2] regions of lower 
filling factor also represent regions of weaker compressive, ten- 
sile, and shear strengths, whereas regions of higher filling factor 
are stronger Therefore, as an analogue to the Weibull distribu- 
tion we propose an initial distribution of the filling factor given 
by the Gaussian function 



1 



1 

■ exp I 



(17) 



(14) where n(<p) is the number density for a filling factor 



IS 

the width of the distribution function, and (pj^ is the median fill- 
ing factor. The pressure limits of compressive, tensile, and shear 
strengths 2(0), T{<p), and Y{(p), respectively, mark the transition 
from the elastic to the plastic hydrostatic regime. Therefore, they 
can be regarded as analogues to the activation threshold ead- 
Hence, we can associate with each <p analogues to the activa- 
tion threshold, which accounts for the first criterion of a damage 
model (see end of Sect. 13.1b that is based on the distribution of 
defects. 

Damaged regions in this simple inhomogeneity scheme are 
therefore represented by areas of low filling factor. Under con- 
stant tension, becomes increasingly smaller with time in these 
regions and, owing to the 0-dependence of either the tensile 
strength T{(p) or shear strength Y{(p), the strength decreases in 
these regions as in Eq. [14] with increasing D. This fulfils the sec- 
ond criterion: the stress reduction. 

Instead of artificially introducing a damage parameter D and 
an activation threshold fact, the filling factor (p naturally takes 
over this twofold role. On the one hand, it represents a quan- 
tity which determines the activation thresholds for the plastic 
regime, and on the other, it also represents a damage parameter 
that reduces, which decreases the strength quantities. Therefore, 
in the inhomogeneity approach, the time evolution of the ap- 
proximated damage is given by the time evolution of or equiv- 
alently the density. The propagation speed of the defects is, in- 
trinsically, the sound speed. 

We introduce a typical size of an inhomogeneous clump 
to comply with the third criterion. In terms of this parameter, 
an absolute length scale is introduced into the inhomogeneity 
approach. 

A great advantage of our approach is that the material param- 
eters (pa-, <Pi.i, and Rc can be deter mined in laborator y measure- 
ments, e.g. by X-ray tomography (iGiittler et al.ll2009l) . However, 
no systematic study has been performed yet. A single measure- 
ment for a random ballistic deposition dust sample, which has 
the lowest filling factor {(p - 0.15) and the highest degree of ho- 
moge neity, has found a s tandard deviation of (p,r ~ 0.013 (Figure 
10 in lGiittleretai]l2009h . There exist no data on either the typ- 
ical clump size or clump size distribution. Therefore, as a first 
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Fig. 1. Initial filling factor distributions of the target and projec- 
tile as binned (top) and cumulative diagram (bottom). The num- 
ber density n{(p) of volumes with filling factor (p is plotted over 
(f>. A larger degree of inhomogeneity is characterised by the in- 
creasing standard deviation (p,j- of the Gaussian (Eq. [TTl i around 
the initial filling factor (p; of the homogeneous aggregate. Here, 
we choose 0i = 0.35. The values of 0o-, can be found in Table[T] 

approximation, in Sect. |4]we study the effect of (per on the out- 
come of pre-planetesimal collisions for a fixed clump size R^. 
The influence of the latter will be studied in future work. 



3.3. Implementation issues 

The inhomogeneity is imposed on the simulated aggregate as an 
initial condition. Initially, we generate a homogeneous aggre- 
gate by assigning a constant initial filling factor <pi. As a sec- 
ond step, <pi is modified according to a Gaussian distribution 
with a standard deviation <pcr, which we will simply call either 
the "inhomogeneity" or "degree of inhomogeneity" in this pa- 
per. For the distribution, we randomly select a particle a and 
set its density to a new (p" following the Gaussian. The same 
(pi" is assigned to all particles within a radius of typical clump, 
which is fixed to = 9.75 mm, the equivalent of one smooth- 
ing length li. Together with a spatial distance of 2.6 mm between 
the SPH particles, one clump is resolved b y ~ 290 particles, 
which is well above the required resolution dGeretshauser et all 
I2010h . Since there are no empirical data for the typical clump 
size yet, as a first approximation we assume that the clump size is 
constant. With the chosen size, the simulated spherical decime- 



tre bodies can be regarded as being assembled from roughly 
spherical centimetre aggregates in sticking processes with differ- 
ent velocities similar to those described by Guttler et al. (201^ 
and lGeretshauser et al.l (1201 ih . However, the effect of the typical 
clump size and resolution on the resulting fragment distribution 
will be studied in future work. 

It is very likely that dust aggregates are compacted durin g 
their evolution in the protoplanetary disc dZsom et al 1 1201 01) . 
in particular if the aggregates are assembled by smaller units 
in sticking processes. For this reason, we do not consider very 
fluffy aggregates with (pi - 0.15 as produced by random ballis- 
tic deposition. We instead simulate spheres whose filling fac- 
tor lies between the minimum (0min - 0.12) and maximum 
(f^max - 0.58) filling factor of the compressive strength relation 
(Eq.|9]l,i.e. <Pi = 0.35. 

The mass of all individual SPH particles, which is a numer- 
ical parameter in the SPH equations, takes a constant value in 
the simulation. This value was computed to be consistent with 
the initial homogeneous aggregate (here (pi - 0.35). Therefore, 
the modification of the density introduces a slight inconsistency 
in the SPH particle distribution. This means that the density 
value of an arbitrary particle differs slightly from the value com- 
puted from the particle number density using the standard SPH 
sum approach. Potentially, this could introduce pressure fluctu- 
ations causing spurious particle motion. However, in our algo- 
rithm the density is determined by solving the continuity equa- 
tion and not by evalu ating the number density of the particles 
(e.g. lMonaghanll2005h . As a consequence this inconsistency is 
marginal. Moreover, the time evolution of a inhomogeneous dust 
aggregate at rest was simulated and no spurious particle mo- 
tion has been detected. In addition, the aggregate displayed no 
signs of instabilities within simulation times much longer than 
the collision timescale. Consequent ly, the presented approach 
is assumed to yield stable aggregates. Evaluating the generated 
SPH particle distribution yields the resulting filling factor dis- 
tribution: the number density n{<p) of the volume fraction with 
filling factor (p is shown in Fig. [T] as a binned and cumulative 
plot for various (p,j-. 

The inhomogeneous clumps with different filling factors are 
clearly visible in Fig.|2] where the density structure of the target 
is colour coded for different <pa-. The packages are picked ran- 
domly but the random pattern is the same for all aggregates to 
ensure that the results are comparable. Owing to the increasing 
(pa-, the regions with (p + (pi gain higher maximum and lower 
minimum filling factors. Consequently, the inhomogeneity pat- 
tern for (pa-i and (pcri (see Tab. [TJ is barely visible in Fig. |2] (first 
row, second and third). In contrast, for (peri the difference be- 
tween the extreme values of the filling factor create a clearly 
visible pattern with high contrast. 

The standard deviation of the Gaussian filling factor curve 
^cr is supplied to the distribution routine as an input parame- 
ter. The range of (p^ begins slightly below the empirical value 
for a homogeneous random ballistic deposition aggregate {(per - 
0.013). There is no empirical data on the inhomogeneity of ag- 
gregates resulting from subsequent collisions. However, owing 
to the collisional history of an aggregate with macroscopic pre- 
planetesimals it is very likely that realistic aggregates have a 
higher degree of inhomogeneity. Therefore, we increase in 
steps of 0.01 up to (per - 0.050. The real standard deviation com- 
puted from the full width at half maximum (FWHM) is smaller 
than the (pa- supplied to the SPH particle initialisation algorithm 
(algorithmic value). The FWHM values are listed in Tab. [1] For 
the sake of simplicity, we use the algorithmic value this paper. 
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Fig. 2. Interior and exterior view of the targets for different inhomogeneities with the filling factor colour coded. The inhomogeneity 
increases from the homogeneous target (a) to the most inhomogeneous target (f). In particular, the standard deviations are (p,y - 
(a), ^a- - 0.01 (b), 0O- = 0.02 (c), <p,j - 0.03 (d), 0o- = 0.04 (e), and <p,j - 0.05 (f). All targets have a similar pattern that originates 
from considering interacting particle passages in the implementation. With increasing (p^ the maximum and minimum filling factor 
of the different spots increase and decrease, respectively. The projectiles have a similar appearance. The result of collisions among 
these aggregates with 10 ms"' are depicted in Fig. [3] 



Symbol 


algorithmic value 


FWHM value 




0.010 


0.0085 


<l>cr2 


0.020 


0.018 


<^tr3 


0.030 


0.027 


<l>a4 


0.040 


0.034 


(pa-5 


0.050 


0.042 



Table 1. Selection of the standard deviations (pu- as they are used 
in the inhomogeneity algorithm. They are compared with the 
achieved value deduced from the FWHM of resulting filling fac- 
tor distribution in Fig. [T] 



4. Simulation results 

We now perform test simulations using our inhomogeneity ap- 
proach. The initial setup is the same as in most simulations of 
iGeretshauser et all (|20I ih . It consists of two colliding spheres 
with <pi = 0.35. The radii of the target and projectile are rt - 
10cm and rp = 6cm, and they consist of 238,238 and 51,477 
SPH particles, respectively, with masses of 1.23 X 10"^ kg each. 
The SPH particles are set on a cubic grid separated by a distance 
of 2.6mm. We consider two impact velocities of vq = 10 ms"' 
and I'o = 12.5 ms"'. For homogeneous aggregates, the smaller 
velocity leads to the sticking of the projectile to the target (see 
Fig. El a). In contrast, the higher velocity results in fragmenta- 
tion (see also Geretshauser et al. 201 1). For these two collision 
velocities, the inhomogeneity of both the target and the projectile 
is defined by the standard deviation of the Gaussian (pa-, which 



is varied from to 0.05 yielding the distributions depicted in 
Fig.[T] The resulting targets are shown in Fig.|2] The projectiles 
display a similar density pattern. 

For vo = 10ms"', the final fragment distribution is shown in 
Fig. [3] from the impact direction. For homogeneous aggregates 
(a) and the small standard deviation (0o- - 0.01, b), the target 
remains intact and the target and projectile form one large ag- 
gregate. For (per - 0.01, some small fragments are visible. For 
(pa ^ 0.02, the target fragments and with greater inhomogeneity 
the largest fragments of the distribution decrease in size. 

The first impression from our visual inspection of the 
collision products is confirmed by the more detailed analy- 
sis. For this purpose, we utilise the four-population model 
(IGeretshauser et all |20I ih . According to this classification 
scheme of collision outcome, four classes of fragments are dis- 
tinguished according to their mass: the largest fragment (in- 
dex "1") and the second largest fragment (index "2") are char- 
acterised by single values of mass m, kinetic energy E, and 
other parameters. To account for the resolution limit, the sub- 
resolution population (index "sr") is introduced, which contains 
all fragments consisting of a single SPH particle (mass mspe)- 
This class is mainly described by global quantities. Fragments 
with masses mspH < »ifrag < m2 form a class whose cumula- 
tive mass distribution is described by a power-law (index "pw"). 
This class is characterised by global quantities and distribution 
functions. 

The evolution of the masses of each population with (pa- 
is shown in Fig. |4] as a fraction of the total mass. The upper 
and the lower plots display the results for vq = 10 ms"' and 
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Fig. 3. Outcome of a collision between a target and projectile with the same (^o- for different inhomogeneities. In all cases, the 
target and projectile radii are ri - 10 cm and rp = 6 cm, respectively, and the collision velocity is vq = 10ms"'. Analogous to 
the initial targets in Fig.|2] the standard deviations for the Gaussian are 0o- = (a), (pa- - 0.01 (b), (p^ - 0.02 (c), (p„- - 0.03 (d), 
(p,r = 0.04 (e), and 0o- = 0.05 (f). The collision outcome is shown in the impact direction. In the homogeneous case (a) and for small 
inhomogeneities (b), the target stays intact and forms one massive object with the projectile. For 0o- > 0.02, the target fragments 
(c-d). The fragment sizes decrease with increasing (per and at the same time the number of fragments increases. 




Fig. 4. The total masses of the four populations with inhomo- 
geneities measured by <pcr- The collision velocities are vo = 
10ms"' (top) and vo = 12.5 ms"' (bottom). The masses of the 
largest fragment mi, of the second largest fragment 1112, of the 
power-law population mp„, and of the sub-resolution popula- 
tion are stacked and normalised by the total mass of the 
system mtot. In both cases, nii and m2 decrease with increas- 
ing (pa-. Conversely, mp„ and m^^ increase. For the homogeneous 
case {(pa- = 0), sticking is found for the low velocity collision, 
whereas for vo = 12.5 ms"' 1112 + and m^^ + 0, which indi- 
cates fragmentation. 



Vo = 12.5 ms" , respectively. The exact values can be found in 
Tab. 12] For the lower collision velocity and ^a ^ 0.01, nearly all 
mass is stored in the largest fragment m\. For (pa > 0.01, a sec- 
ond largest fragment and a power-law population appears. With 
increasing inhomogeneity mi rapidly decreases, while the mass 
of the second largest fragment m2 only slightly decreases. The 
masses mi and 1112 become comparable in the end. For high in- 
homogeneity values, most of the mass is stored in the power-law 
population nip^, as we discuss below. However, for high inhomo- 
geneities mpw seems to remain fairly constant, while the mass 
of the sub-resolution population msr increases. For the higher 
collision value (bottom plot of Fig. HJ, the evolution is similar 
but starts with mpw 0. The mass mi constantly decreases, 
but not as rapidly as in the low velocity case. Moreover, m2 
at first increases, then slightly decreases until the largest and 
second largest fragment become comparable. The mass of the 
power-law population rapidly increases for 0.01 < (pa < 0.04. 
For larger (pa, the power-law population only slightly increases, 
whereas msr increases at a higher rate. The mass evolution with 
increasing inhomogeneity can be summarised as follows: (1) in- 
homogeneity leads to fragmentation at collision velocities, at 
which homogeneous aggregates do not fragment, i.e. inhomo- 
geneous aggregates are weaker (2) For the two investigated ve- 
locities a larger (pa leads to a decrease in the mass of the largest 
and second largest fragments. (3) A higher degree of inhomo- 
geneity causes an increase in both the masses of the power-law 
and sub-resolution population. 

We turn to the analysis of the distribution of the residual total 
energy after the collision, which is shown in Fig.|5]as a fraction 
of the initial kinetic energy. The residual total energy is the sum 
of the translational energy and the energy stored in the inter- 
nal degrees of freedom (e.g. rotation, libration) of the respec- 
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Fig. 5. The energy of each population for increasing inhomo- 
geneity parameter 0o- and collision velocity vq = 10 ms"' (top) 
and I'o = 12.5 ms"' (bottom). The energy of each population 
is the sum of translational, rotational and librational energy and 
normalised by the initial kinetic energy Ei„it . The two latter en- 
ergy contributions are negligible for the head-on collisions pre- 
sented here. Owing to the strong correlation with mass, this di- 
agram follows a similar trend to Fig. |4] The energy contribu- 
tions are divided up into £1, E2, Ep„, and £51 for the largest 
and second largest fragment, the power-law and sub-resolution 
population, respectively. In both velocity cases, both Ei and £2 
decrease with increasing (pg-. In contrast, Ep„ and £51 increase. 
In particular, the high fraction of kinetic energy stored in the 
sub-resolution population is remarkable. 



live populations. The tabulated values can be found in Tab. |2] 
In Fig.|5] the low velocity case with vq = 10 ms"' is shown in 
the upper plot and the high velocity case with vo = 12.5 ms"' 
in the lower plot. For the former, a large amount of energy £1 
is stored in the largest fragment, but a considerable amount £si 
is also stored in the sub-resolution population. As £1 decreases, 
the contribution of the second largest fragment at first increases, 
then decreases slightly with larger <pij. As for the mass, the en- 
ergy contribution of the power-law population £p„ increases 
with a higher degree of inhomogeneity. Remarkably, the total en- 
ergy fraction of these three populations remains almost constant. 
In contrast, the energy contribution of the sub-resolution popu- 
lation £si increases more strongly than its mass contribution in 
Fig. |4] The overall residual energy increases with a higher de- 
gree of inhomogeneity and this energy excess is stored mostly in 
the sub-resolution population. For the high velocity case (bottom 
plot of Fig.jSj, this is even more evident: £1, £2, and £pw behave 
similar to the low velocity case but with £2^0 and £pw for 
= 0. For a higher degree of inhomogeneity, less energy can 
again be dissipated and the total residual energy of the largest 
three populations remains nearly constant for increasing (pa-- As 
a consequence, the residual energy excess is mainly stored in 
the sub-resolution population for large values of <pa-. The conclu- 
sions of our energy analysis are; (1) with increasing impact ve- 
locity less energy can be dissipated by the sy stem and the resid- 
ual en ergy increases as previously found by iGeretshauser et all 
(1201 Ih . (2) With increasing inhomogeneity, the eff'ectiveness of 



Fig. 6. Cumulative mass of the fragment masses over their aver- 
age filling factor (p. The mass is normalised by the total fragment 
mass ffitot- The collision velocities are vq - 10 ms"' (top) and 
Vo = 12.4 ms"' (bottom). Compared to the initial distribution 
of filling factors (Fig. [TJ, the curves are shifted significantly to- 
wards higher filling factors. This reflects the compression taking 
place during the collision. 



the energy dissipation decreases for the two investigated veloci- 
ties. (3) The energy contribution of all populations show a sim- 
ilar dependence on to their mass contributions: the fractions 
of energy stored in the largest and second largest fragment fol- 
low a decreasing trend with increasing inhomogeneity. The con- 
tribution of power-law and sub-resolution populations increase. 
(4) While the overall energy contribution of the largest fragment 
populations remains nearly constant with 0o-, the fraction stored 
in the sub-resolution population drastically increases, such that 
the energy that cannot be dissipated owing to the increasing in- 
homogeneity is stored in rapidly moving single SPH particles. 

In comparison to the initial filling factor distribution of 
Fig.[T] the final filling factor distribution of the cumulated mass 
is shifted towards higher filling factors. Figure |6]reveals that in 
most cases more than 85 % of the fragment mass is stored in fill- 
ing factors > 42 % (except for vq - 12.5 ms"' and <pa- - 0.050). 
With the same exception, the distributions look quite similar in 
both velocity cases. There is a slight tendency for there to be 
more mass at lower filling factors in the collision with \\) = 
12.5 ms"'. From Fig. |6l it can also be seen that the broader the 
initial filling factor distribution, the broader the final filling fac- 
tor distribution. For the vq = 10.0 ms"' case, the two most ho- 
mogeneous aggregates hardly produce any fragments that cause 
a jump in the curves without any intermediate steps. The dis- 
tribution for the highest collision velocity (vq - 12.5 ms"') and 
the highest degree of inhomogeneity ((per = 0.050) is particularly 
interesting since ~ 30 % of the fragment mass possess filling fac- 
tors < 0.35, which is a larger fraction than for the initial aggre- 
gate. This implies that rarefaction occurs during the collision. On 
the microscopic scale, this is equivalent to monomer re-ordering 
and stretching of monomer chains. We note that a tiny fraction 
of particles acquire overdensities {(f> > 1) in the collision. In the 
most violent case, this applies to single SPH particles, which in 
total represent 0.03 % of the total mass of the system, which is 
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Fig. 7. Cumulative mass distribution of a power-law population 
for different inhomogeneity parameters <pa- and collision veloci- 
ties I'o - 10 ms"' (top) vo - 12.5 ms"' (bottom). In both veloc- 
ity cases both the mass of the largest member of the power-law 
population and the power-law slope decrease. This indicates the 
fragmentation to smaller aggregates for increasing (pi^. For equal 
(pa-, the slopes are shallower for the higher velocity. The power- 
law fit parameters can be found in Tab. [3] 



negligible. These overdensities indicate that the monomers are 
breaking up on the microscopic scale. 

We have so far considered all members of the four- 
population model. We now constrain our analysis to the power- 
law population. Its mass distribution is shown in the cumulative 
plot of Fig. I?] again for Vo - 10ms"' (top) and vo = 12.5 ms"' 
(bottom). The figure shows the ratio of the cumulative mass mcum 
to the fragment mass mftag, which are both normalised by the 
total mass of the power-law population nip^ giv en in Tab. |2] 
The c umulative distribution can be accurately (e.g. lGiittler et al.l 
120 Id) described by a power law 



n('Wfrag) = J «('«) m dm = ^ 



»Ifrag 



(18) 



where n(m) dm is the number density of fragments of mass m, 
which is normalised by the total mass of the power-law pop- 
ulation ;7ip„. The quantity jUp„ is the normalised mass of the 
most massive member of the power-law population and k is the 
power-law index or fragmentation parameter. We fitted the cu- 
mulative mass distributions of Fig. Q for all fragment masses 
with Eq. [18] Finding an explanation for the kink in the cumu- 
lative mass curves is the subject of ongoing research. This devi- 
ation from a po wer law might be either a physical or a resolution 
effect (see also iGeretshauser et al.ll20in) . The results are listed 
in Tab. [3] Figure [8] shows that the power-law index k decreases 
slightly as the degree of inhomogeneity increases for both colli- 
sion velocities, which leads to increasingly shallower slopes in 
Fig. Blowing to the production of smaller fragments with increas- 
ing In addition, the k values for the higher impact velocity 
(vo = 12.5 ms"') are systematically below the k values of the 
lower velocity (v'o - 10.0 ms"'), except for small values of (per 
for which there is only a small number of fragments (see Tab.[2]i 
and hence the quality of the statistics is insufficient. This find- 
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Fig. 8. Fragmentation parameter k of the power-law population. 
For increasing inhomogeneity, indicated by the standard devia- 
tion (pa, the fragmentation parameter decreases. This indicates 
shallower slopes of the power-law mass distribution which is 
interpreted as an increasing fraction of smaller fragments. For 
(pa- ^ 0.02, the higher impact velocity results in smaller values 

of K. 







K 







10.0 


0.44 ± 0.027 


0.10 ± 2.4 X 10"' 


0.01 


10.0 


0.88 ± 0.035 


0.10 ±2.3 X 10"' 


0.02 


10.0 


0.88 ± 0.046 


0.36 ±0.012 


0.03 


10.0 


0.59 ±0.011 


0.24 ± 6.9 X 10"' 


0.04 


10.0 


0.60 + 6.4 X 10"^ 


0.11 ± 1.7 X 10"' 


0.05 


10.0 


0.38 ± 5.4 X 10"' 


40. 10 ± 4.0 X 10"' 





12.5 


3.88 ±0.254 


0.20 ± 1.4 X 10"' 


0.01 


12.5 


0.92 ± 0.030 


0.52 ± 1.0 X 10"- 


0.02 


12.5 


0.54 ±0.017 


0.18 ± 1.0 X 10"^ 


0.03 


12.5 


0.46 ± 0.006 


0.10 ± 3.0 X 10"' 


0.04 


12.5 


0.41 ±0.004 


0.06 ± 1.3 X 10"' 


0.05 


12.5 


0.23 ± 0.003 


0.02 ± 1.0 X 10"-' 



Table 3. Fit values of the power-law population. The quantities k 
and yUpw denote, respectively, the slope of the power-law fit and 
the mass of the most massive member, which is normalised by 
the total mass of the power-law population wipw. The standard 



deviation of the inhomogeneity Gaussian is denoted by 
the impact velocity by vq. 



and 



ing indicates that a larger fraction of less massive fragments are 
produced at higher impact velocities, which is expected. 

As shown in Fig. |9] the mass of the largest member of the 
power-law population, given by //pw, at first increases for larger 
values of (pa-. It has a maximum at (pa ~ 0.02 for vq = 10 ms"' 
and at (pa- ~ 0.01 for vq = 12.5 ms"', and then decreases in 
an exponential fashion. When comparing both velocity cases 
(see Tab. |2]i, we find that more fragments are produced in the 
high v elocity case. This was also observed bvlGeretshaus erTt al.l 
(|20I ih . The influence of the inhomogeneity on the power-law 
population can be summarised as follows: (1) the overall mass 
of the power-law distribution is increased as discussed above; 
(2) a higher degree of inhomogeneity leads to the production of 
smaller fragments; (3) the largest member of the power-law pop- 
ulation reaches its highest mass at small (pa-, after which its mass 
decreases with increasing (pa-. 
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Table 2. Simulation results presented according to the four-population classification. The subscripts "1", "2", "pw", "sr", and "tot" 
denote the largest and second largest fragment, the power-law, sub-resolution, and total fragment population, respectively. The 
quantities E, m, and are the total energy, the mass, and the number of fragments, respectively. The standard deviation for the 
inhomogeneity Gaussian distribution is given by (pa- and the impact velocity by vq. 
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Fig. 9. Fragmentation strength of the power-law population. The 
quantity /ip„ represents the normalised mass of the most massive 
member of the power-law population and serves as an indica- 
tor of the fragmentation strength. For increasing inhomogeneity, 
represented by larger values of (pa-, fJ-pw at first increases and then 
decreases again. The increase may be caused by our low quality 
statistics owing to the small number of fragments. 



5. Discussion and outlook 

We have presented the first simulations of inhomogeneous 
porous pre-planetesimals, where the inhomogeneity of SiOi dust 
aggregates is explicitly resolved. This inhomogeneity can 
in principal also be determined in laboratory experiments 
dGiittler et al.ll2009l) . The approa ch is based on the conce pt that 
according to our porosity model (iGeretshauser et al.l2010l) inho- 
mogeneities in filling factor cause fluctuations in the compres- 
sive, shear, and tensile strength in the aggregate. These fluc- 
tuations can be regarded as flaws in the mater ial. In contrast 
to da mage models designed for brittle material dGradv & KippI 
Il980l) . we did not explicitly evolve the propagation of these 
flaws. Instead, the defects in the dust material, which behaves 
more like a fluid, were approximated by the time evolution of the 
filling factor or - equivalently - the density. The inhomogeneity 
of an aggregate was assumed in the initial SPH particle distribu- 



tion by using a Gaussian distribution to model the filling factor. 
The inhomogeneity is characterised by the standard deviation (pu 
of the Gaussian and the typical radius of a clump. 

Using this inhomogeneity approach, we have performed 
test simulations of collisions between dust aggregates of in- 
termediate porosity with fixed typical clump size. Two colli- 
sion velocities have been chosen: one below and one above the 
fragmentation threshold for homogeneous aggregates. The re- 
sults have been analysed using the four-population model of 
IGeretshauser et all ( 1201 Ih . For the lower collision velocity, we 
have shown that inhomogeneity leads to fragmentation. For 
both velocities, the masses of largest and second largest frag- 
ment are lower for a higher degree of inhomogeneity whereas 
the masses of the power-law and sub-resolution population are 
higher Focussing on the power-law population, the number of 
fragments and the fraction of small fragments increase with in- 
creasing These findings show that the inhomogeneity ap- 
proach has passed the first test in explicitly resolving the het- 
erogeneity of porous objects. In a future study, we will address 
the influence of the typical clump size on the outcome of pre- 
planetesimal collisions. 

Our findings indicate that inhomogeneous dust aggregates 
are weaker than their homogeneous equivalents. A slight inho- 
mogeneity is sufficient to cause a catastrophic disruption instead 
of growth as the result of a dust aggregate collision. Therefore, 
inhomogeneity might explain the higher velocity thresholds for 
fragmentation in sim ulations than the a lower v alue found in lab- 
oratory experiments (iGeretshauser et al]l201 ih . 

Furthermore, macroscopic dust aggregates in protoplanetary 
discs are produced by the subsequent impacts of smaller aggre- 
gates at different impact velocities, i.e. pre-planetesimals have a 
collision history, thus are very likely to be inhomogeneous. This 
has a negative effect on planetesimal formation by subsequent 
collisions because pre-planetesimals might become weaker as 
the number of collisions increases. As a result, their threshold 
velocity for catastrophic disruption might be smaller. However, 
in a future study we will investigate the influence of the clump 
size of the stability of the pre-planetesimals. If the clumps are 
small relative to the pre-planetesimal size, this might increase 
the stability of the aggregates. 

With the model presented in this paper, we are able to include 
these features. Further studies might be carried out that investi- 
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gate the inhomogeneity created by subsequent multiple impacts. 
The result could be classified according to the inhomogeneity 
model that we have presented here. 

To date it has been found that the simulations of collisions 
between homogeneous agg regates carried o ut by means of the 
porosity model by Geretsha user et al.l (1201 Ol) depend only on the 
filling factor and produce the same fragment distribution for an 
identical set of input parameters. By randomly assigning an in- 
homogeneity pattern it is now possible to profoundly investigate 
the statistics of a fragment distribution. Statistical fluctuations in 
the quantities of the four-population model can now be estimated 
in particular for simulations with a small number of fragments. 

In high velocity grazing collisions, the target begins to rotate. 
For highly porous aggregates, which have a low tensile strength, 
high spinning rates tend to lead to fragmentation of the target. 
With increasing inhomogeneity, this might also be true for ag- 
gregates with medium and low porosity. A quantitative investi- 
gation of this effect can be carried out by means of the presented 
approach. 

The filling factor distribution of dust aggregates can be de- 
termined in the labo ratory by X-ray tomography measurements 
dGiittler et al.ll2009l) . These empirical data can be directly imple- 
mented into our inhomogeneity mod el. In additio n the success- 
ful material calibration (Geretshause r et aDl2010h . the obtained 
inhomogeneity parameters can be used to improve the realistic 
simulation of porous dust aggregates. 
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^tot 
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h 
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second irreducible invariant 
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Ko 


bulk modulus for the uncompressed dust sample 


m 
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r masses of largest and second largest fragment. 


power-law and sub-resolution populations 


^cum 


cumulated mass 
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projectile mass 
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target mass 
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exponent in the Weibull distribution 
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number density of volumes with filling factor ^ 
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strain activation threshold 




strain tensor 
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p 
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matrix density 




stress tensor 


<P 
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0max5 0min 
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initial filling factor 




inhomogeneity median 
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compressive strength 
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